We study the systematics effects induced by gain mismatches for 1 yr Litebird scanning strategy and a fake focal plane with just 2 detector pairs, (total 4 detectors). Also, I 've run a couple of tests simulating South Pole data with Atmospheric sims ( -> Reijo ).
scramblegains routine ) only in one of the two detectors ( namely $detB$) as a random Gaussian fluctations around 0 and with a 5% width during a fraction of precession period ($5, 10, 96$ min). We refer to the TODs with gaindrifts injected as pre-correction, whereas the ones corrected with the methodology described above as post-correction.
Simulate 1 yr of observations with LB scanning strategy parameters, α=45, β=50, spinrate=0.1rpm, precperiod=96.5min Focal plane encodes 4 detectors at 140 GHz , FWHM= 36 arcmin, NET= 38.6 μK√s, sample rate= 40 Hz and not spinning HWP.
we firstly want to quantify how the gain correction procedure is affected by the template chosen to perform the regression . We thus consider signal only simulations encoding as a template for the regression either the Orbital Dipole or the CMB Temperature.
We assess the T->P leakage both at the map and i power spectrum level.
In the case of signal only sims, there is no need to destripe. We therefore looked at binned maps.
from IPython.display import HTML
HTML('''<script>
code_show=true;
function code_toggle() {
if (code_show){
$('div.input').hide();
} else {
$('div.input').show();
}
code_show = !code_show
}
$( document ).ready(code_toggle);
</script>
<form action="javascript:code_toggle()"><input type="submit" value="Click here to toggle on/off the raw code."></form>''')
import healpy as hp
import pylab as pl
inputclfile ='/Users/peppe/work/satellite_sims/sigonly/inputcl.fits'
try :
cl = hp.read_cl(filename=inputclfile)
except IOError:
cl=hp.anafast(cmbmaps, lmax=lmax)
hp.write_cl(cl=cl, filename='inputclfile' )
hitmap=hp.read_map('/Users/peppe/work/satellite_sims/madam_hmap.fits', verbose=False)
hp.mollview(pl.log10( hitmap), min=3.5, max=4, unit='Log (hitcounts)' ,title='')
cl = pl.loadtxt('/Users/peppe/work/satellite_sims/Class_planck2015_r0_cl_lensed.dat')
lmax=3*hp.get_nside(cmbmaps[0])
l=pl.arange(2, lmax+1)
Nl =lambda w: (pl.pi /180. /60. * w)**2
wp=7.6 # uK arcmin
Nl0= Nl(wp)
dipolemap= hp.read_map(field=range(3) ,
filename='/Users/peppe/work/satellite_sims/sigonly/dipole/madam_bmap.fits',
verbose=False)
sigonlymap = hp.read_map(field=range(3) ,
filename='/Users/peppe/work/satellite_sims/sigonly/dipole_CMB_T/madam_bmap.fits',
verbose=False)
fig=pl.figure(figsize=(15,7))
hp.mollview(sigonlymap [0], unit='K ', title='Dipole T ', sub=221 , fig=fig )
cmbT =sigonlymap [0] - dipolemap[0]
hp.mollview(cmbT , unit='K ', title='CMB T ', sub=222 , fig=fig )
hp.mollview(sigonlymap [1] , unit='K ', title='Q ', sub=223 , fig=fig )
hp.mollview(sigonlymap [2], unit='K ', title='U ', sub=224 , fig=fig )
driftmap= hp.read_map(field=range(3) ,
filename='/Users/peppe/work/satellite_sims/sigonly/uncorr/madam_bmap.fits',
verbose=False)
fig=pl.figure(figsize=(15,8))
hp.mollview( driftmap [0], unit='K ', title=r'$ T$ ', sub=131 , fig=fig )
hp.mollview(driftmap [1], unit='K ', title=r'$ Q$', sub=132 , fig=fig )
hp.mollview(driftmap [2], unit='K ', title=r'$U$', sub=133 , fig=fig )
cldriftdip=hp.anafast(driftmap , lmax=lmax)
Caption: We get similar non-null signal on Q, and U maps: that's the expected T-> P leakage from gain drifts. In T map the leakage is $\sim 2$ orders of magnitude above the input signal (CMB T + Dipole ). A ring-like pattern is observed in Q and U maps and seems to fall exactly where the dipole changes sign.
We now estimate the gains by regressing against the Dipole template.
driftcorrmap= hp.read_map(field=range(3) ,
filename='/Users/peppe/work/satellite_sims/sigonly/corr/madam_bmap.fits',
verbose=False)
fig=pl.figure(figsize=(15,8))
hp.mollview((driftcorrmap [0] ) ,unit='K', title=r'$ T $ ', sub=131 ,
fig=fig )
hp.mollview(driftcorrmap [1] , unit='K',title=r'$ Q $', sub=132 , fig=fig )
hp.mollview(driftcorrmap [2] ,unit='K', title=r'$ U $', sub=133 , fig=fig )
cldriftcorrdip=hp.anafast(driftcorrmap , lmax=lmax)
Caption: Notice that Q and U maps are a factor ~10 times smaller than the pre-corr. ones (see above ). They present both a dipole-like structures and a couple of scanning strategy features (the two horizontal stripes along the Ecliptic plane and in the Ecliptic poles) . The latter are the regions where pixels are observed $\sim 3 $ times more wrt the rest of the sky, as shown in the hitmap above.
fig=pl.figure(figsize=(16, 5))
pl.subplot(121)
pl.title('EE')
pl.ylabel(r'$C_{\ell} [ \mu K^2 ]$', fontsize=16)
pl.loglog(l, cl[ 2:lmax+1,2] ,'k', label='CMB')
pl.loglog(l, cldriftdip[1][2:] *1e12, label='pre-correction')
pl.loglog(l, cldriftcorrdip[1][2:] *1e12, label='post-correction')
pl.legend(loc='upper right')
pl.xlim([2,lmax])
pl.subplot(122)
pl.title('BB')
pl.loglog(l, cl[ 2:lmax+1,3] ,'k')
pl.loglog(l, cldriftdip[2][2:] *1e12)
pl.loglog(l, cldriftcorrdip[2][2:] *1e12)
pl.loglog(l, pl.ones_like(l)*Nl0,'r', alpha=.5, linewidth=3)
pl.ylim([1e-8,1e0])
foo=pl.xlim([2,lmax])
Caption: Solid blue (orange) are power spectra of pre(post)-gain-correction residual maps. Since no noise has been coadded to the TODs so far, the plots above represent the level of accuracy our estimator might achieve whith the Dipole as a template. As a reference, we plot the CMB spectra from $\Lambda CDM$. We also include in BB plot the solid red line resembling the Litebird whitenoise level of $7.6 \mu $K arcmin, after FG subtraction.
Pros: By comparing the blue with the orange solid lines should highlights the effects of gain correction. As expected they are small for the temperature maps, but in terms of EE and BB power the correction is very effective: it may correct the T->P leakage as much as 3 orders of magnitude in the $\sim$ degree scales.
Cons: the correction doesn't work at $\ell<10$ scales. The reason is related to the presence of Dipole: in facts, it introduces $\ell \geq 2 $ bias in the CMB that needs to be removed by pairing a dipole-fitting code with comp.sep. algorithms (see Litebird CDR, sect. 4.6.1).
We keep the same simulation setup as before, we just changed the input maps encoding now only CMB T (assuming Dipole was correctly subtracted).
sigonlymap = hp.read_map(field=range(3) ,
filename='/Users/peppe/work/satellite_sims/sigonly/CMB_T/madam_bmap.fits',
verbose=False)
fig=pl.figure(figsize=(12,8))
hp.mollview(sigonlymap[0] , unit='K ', title=' T ', sub=131 , fig=fig )
hp.mollview(sigonlymap [1] , unit='K ', title='Q ', sub=132 , fig=fig )
hp.mollview(sigonlymap [2], unit='K ', title='U ', sub=133, fig=fig )
driftmap= hp.read_map(field=range(3) ,
filename='/Users/peppe/work/satellite_sims/sigonly/uncorr/wodipole/madam_bmap.fits',
verbose=False)
fig=pl.figure(figsize=(15,8))
hp.mollview( driftmap [0] , unit='K ', title=r'$ T$ ', sub=131 , fig=fig )
hp.mollview(driftmap [1], unit='K ', title=r'$ Q$', sub=132 , fig=fig, norm='hist' )
hp.mollview(driftmap [2], unit='K ', title=r'$ U$', sub=133 , fig=fig , norm='hist' )
cldriftsT=hp.anafast( driftmap , lmax=lmax)
Caption: The pre-correction maps look as noise fluctuating around $10^{-6} $ K. Note that we've shrunk the colorscale to highlight differences .
We now apply the gain correction and output the maps.
driftcorrmap= hp.read_map(field=range(3) ,
filename='/Users/peppe/work/satellite_sims/sigonly/corr/wodipole/madam_bmap.fits',
verbose=False)
fig=pl.figure(figsize=(15,8))
hp.mollview(driftcorrmap [0] , unit='K ', title=r'$ T$ ', sub=131 , fig=fig )
hp.mollview(driftcorrmap [1] , unit='K ', title=r'$ Q$', sub=132 , fig=fig, norm='hist' )
hp.mollview(driftcorrmap [2], unit='K ', title=r'$ U$', sub=133 , fig=fig, norm='hist' )
cldriftcorrT=hp.anafast(driftcorrmap , lmax=lmax)
Caption: Post-correction residual maps are one order of magnitude smaller than the pre-correction ones, indicating that the methodology as a clear effect at the map level. Also, we achieve 10 times lower residuals by regressing with CMB T maps than in the Dipole case. A T-like pattern is visible both in Q and U maps .
fig=pl.figure(figsize=(16, 5))
pl.subplot(121)
pl.title('EE')
pl.ylabel(r'$C_{\ell} [ \mu K^2 ]$', fontsize=16)
pl.loglog(l, cl[ 2:lmax+1,2] ,'k', label='CMB')
pl.loglog(l, cldriftsT[1][2:] *1e12, label='pre-correction')
pl.loglog(l, cldriftcorrT[1][2:] *1e12, label='post-correction')
pl.legend(loc='upper right')
pl.xlim([2,500])
pl.ylim([1e-8,1e-2])
pl.subplot(122)
pl.title('BB')
pl.loglog(l, cl[ 2:lmax+1,3] ,'k')
pl.loglog(l, cldriftsT[2][2:] *1e12)
pl.loglog(l, cldriftcorrT[2][2:] *1e12)
pl.loglog(l, pl.ones_like(l)*Nl0,'r', alpha=.5, linewidth=3)
pl.ylim([1e-8,1e-2])
foo=pl.xlim([2,500])
Caption: Power spectra pre and post correction in the case of T only sims. Notice that the systematic constamination is lower as compared to the one with Dipole, though even large scale B-modes ($\ell <10 $) of post-correction maps are as much as comparable as the expected level of lensing B-modes at those angular scales. However, we are safely below the LB white noise level after Foreground subtraction (solid red line).
from IPython.display import HTML
HTML('''<script>
code_show=true;
function code_toggle() {
if (code_show){
$('div.input').hide();
} else {
$('div.input').show();
}
code_show = !code_show
}
$( document ).ready(code_toggle);
</script>
<form action="javascript:code_toggle()"><input type="submit" value="Click here to toggle on/off the raw code."></form>''')
We include scalar $E$ -modes polarization in the input maps. Also, TE power spectrum is no more null.
sigonlymap = hp.read_map(field=range(3) ,
filename='/Users/peppe/work/satellite_sims/sigonly/scalar_only/madam_bmap.fits',
verbose=False)
fig=pl.figure(figsize=(14,12))
hp.mollview(sigonlymap[0] , unit='K ', title=' T ', sub=131 , fig=fig )
hp.mollview(sigonlymap [1] , unit='K ', title='Q ', sub=132 , fig=fig )
hp.mollview(sigonlymap [2], unit='K ', title='U ', sub=133 , fig=fig )
clS=hp.anafast(sigonlymap , lmax=lmax)
Caption: TQU maps of scalar anisotropies (no lensing, no tensors).
driftmap= hp.read_map(field=range(3) ,
filename='/Users/peppe/work/satellite_sims/sigonly/uncorr/scalar_only/10min/madam_bmap.fits',
verbose=False)
fig=pl.figure(figsize=(15,8))
hp.mollview( driftmap [0] , unit='K ', title=r'$ T$ ', sub=131 , fig=fig )
hp.mollview(driftmap [1] , unit='K ', title=r'$ Q$', sub=132 , fig=fig )
hp.mollview(driftmap [2] , unit='K ', title=r'$ U$', sub=133 , fig=fig )
cldriftS=hp.anafast(driftmap , lmax=lmax)
Caption: pre-correction TQU maps w/ 10 min gain drifts .
driftcorrmap= hp.read_map(field=range(3) ,
filename='/Users/peppe/work/satellite_sims/sigonly/corr/scalar_only/10min/madam_bmap.fits',
verbose=False)
fig=pl.figure(figsize=(15,8))
hp.mollview(driftcorrmap [0] , unit='K ', title=r'$ T$ ', sub=131 , fig=fig )
hp.mollview(driftcorrmap [1] , unit='K ', title=r'$ Q$', sub=132 , fig=fig )
hp.mollview(driftcorrmap [2] , unit='K ', title=r'$ U$', sub=133 , fig=fig )
cldriftcorrS=hp.anafast(driftcorrmap , lmax=lmax)
Caption: post-correction TQU maps.
We now plot the EE, BB and TE power spectra.
bl=hp.gauss_beam(lmax=lmax, fwhm=pl.deg2rad(36./60))[2:]
fig=pl.figure(figsize=(16, 5))
pl.subplot(121)
pl.title('EE')
pl.ylabel(r'$C_{\ell} [ \mu K^2 ]$', fontsize=16)
pl.loglog(l, cl[ 2:lmax+1,2] ,'k', label='CMB')
pl.loglog(l, cldriftS[1][2:]/bl**2 *1e12, label='pre-correction')
pl.loglog(l, cldriftcorrS[1][2:]/bl**2 *1e12, label='post-correction')
pl.legend(loc='upper right')
pl.xlim([2,600])
pl.ylim([1e-8,1e0])
pl.subplot(122)
pl.title('BB')
pl.loglog(l, cl[ 2:lmax+1,3] ,'k')
pl.loglog(l, cldriftS[2][2:] /bl**2 *1e12)
pl.loglog(l, cldriftcorrS[2][2:] /bl**2 *1e12)
pl.loglog(l, pl.ones_like(l)*Nl0,'r', alpha=.5, linewidth=3)
pl.ylim([1e-8,1e-2])
foo=pl.xlim([2,600])
fig=pl.figure()
pl.title('TE')
ll=l*(l+1)/pl.pi/2.
pl.ylabel(r'$D_{\ell} [ \mu K^2 ]$', fontsize=16)
pl.plot(l, cl[ 2:lmax+1,4] *ll ,'k')
pl.plot(l, cldriftS[3][2:]*ll /bl**2 *1e12)
pl.plot(l, cldriftcorrS[3][2:]*ll /bl**2 *1e12)
pl.ylim([-2e2,2e2])
foo=pl.xlim([2,700])
Caption: Power Spectra of pre and post correction maps encoding scalar anisotropies (temp. and polarization). (left panel) E-mode spectra. For consistency, with CMB theoretical spectrum we applied a 35 arcmin gaussian beam, whose effects are evident at higher multiples (right panel). (Bottom panel) TE power spectra (multiplied by the normalization factor , $D_\ell = \ell (\ell +1) C_\ell/2\pi$. Notice that at $\ell > 600$ the pre-correct. power spectrum deviates from the theoretical one whereas the post-correction doesn't.
from IPython.display import HTML
HTML('''<script>
code_show=true;
function code_toggle() {
if (code_show){
$('div.input').hide();
} else {
$('div.input').show();
}
code_show = !code_show
}
$( document ).ready(code_toggle);
</script>
<form action="javascript:code_toggle()"><input type="submit" value="Click here to toggle on/off the raw code."></form>''')
To better understand the features observed above we 've run TOAST sims on 20 different CMB realizations (no tensor and lensing B-modes).
We also estimate the transfer function defined as :
$$ \mathcal{F}_\ell = \frac{C_\ell^X}{C_\ell^{input}} $$ with $X$ label indicating either pre or post correction power spectra.
meanS,stdS= pl.load('/Users/peppe/work/satellite_sims/sigonly/scalar_only/10min/sigonly_sims.npy')
meanpre,stdpre= pl.load('/Users/peppe/work/satellite_sims/sigonly/scalar_only/10min/sigonly_sims_uncorr.npy')
meanpost, stdpost= pl.load('/Users/peppe/work/satellite_sims/sigonly/scalar_only/10min/sigonly_sims_corr.npy')
bl=hp.gauss_beam(lmax=lmax, fwhm=pl.deg2rad(36./60))[2:]
fig=pl.figure(figsize=(20, 10))
pl.subplot(231)
pl.title('TT')
pl.ylabel(r'$\langle C_{\ell} \rangle_{MC} [ \mu K^2 ]$', fontsize=16)
pl.loglog(l, meanS[0][2:] *1e12/bl**2, label='no drifts', color='k')
#pl.fill_between(l, (meanS[1][2:] -stdS[1][2:])*1e12/bl**2 , (meanS[1][2:] + stdS[1][2:])*1e12/bl**2 , alpha=.4, color='k' )
pl.loglog(l, meanpre[0][2:] *1e12/bl**2, label='pre-correction')
pl.loglog(l, meanpost [0][2:]*1e12/bl**2 , label='post-correction')
pl.xlim([2,600])
pl.ylim([1e-2,1e3])
pl.subplot(232)
pl.title('EE')
#pl.ylabel(r'$\langle C_{\ell} \rangle_{MC} [ \mu K^2 ]$', fontsize=16)
pl.loglog(l, meanS[1][2:] *1e12/bl**2, label='no drifts', color='k')
#pl.fill_between(l, (meanS[1][2:] -stdS[1][2:])*1e12/bl**2 , (meanS[1][2:] + stdS[1][2:])*1e12/bl**2 , alpha=.4, color='k' )
pl.loglog(l, meanpre[1][2:] *1e12/bl**2, label='pre-correction')
pl.loglog(l, meanpost [1][2:]*1e12/bl**2 , label='post-correction')
pl.legend(loc='upper right', fontsize=12)
pl.xlim([2,600])
pl.ylim([1e-5,1e-1])
pl.subplot(233)
pl.title('TE')
ll=l*(l+1)/pl.pi/2.
pl.ylabel(r'$D_{\ell} [ \mu K^2 ]$', fontsize=16)
pl.plot(l, meanS[3][2:]*ll*1e12/bl**2 , label='no drifts', color='k')
#pl.fill_between(l, (meanS[3][2:] -stdS[3][2:] )*ll*1e12/bl**2, (meanS[3][2:] + stdS[3][2:])*ll*1e12/bl**2, alpha=.4, color='k' )
pl.plot(l, meanpre[3][2:]*ll*1e12/bl**2 , label='pre-correction')
pl.plot(l, meanpost [3][2:] *ll*1e12/bl**2, label='post-correction')
pl.ylim([-1.5e2,1.5e2])
foo=pl.xlim([2,700])
pl.subplot(234)
pl.ylabel(r'$\mathcal{F}_{\ell} $', fontsize=20)
pl.semilogx(l,meanpre[0][2:]/meanS[0][2:] ,label='pre-correction')
pl.semilogx(l,meanpost[0][2:]/meanS[0][2:], label='post-correction')
pl.ylim([0,2])
foo=pl.xlim([2,700])
pl.xlabel(r'$\ell$', fontsize=16)
pl.subplot(235)
pl.semilogx(l,meanpre[1][2:]/meanS[1][2:] ,label='pre-correction')
pl.semilogx(l,meanpost[1][2:]/meanS[1][2:], label='post-correction')
pl.ylim([0,2])
foo=pl.xlim([2,700])
pl.xlabel(r'$\ell$', fontsize=16)
pl.legend(loc='best', fontsize=12)
pl.subplot(236)
pl.plot( l,meanpre[3][2:]/meanS[3][2:])
pl.plot(l, meanpost[3][2:]/meanS[3][2:])
pl.xlabel(r'$\ell$', fontsize=16)
pl.ylim([-2,4])
foo=pl.xlim([2,700])
Caption: (top) Mean TT (left) EE (center) and TE (right) power spectra from 20 MC sims. (bottom) Transfer functions for TT (left) EE (center) and TE (right), as defined above. Pre and post correct. transfer functions similarly fluctuates around $1$, as expected. In particular, at $\ell > 400 $, the gain drift systematics inject extra power to both EE and TE (also TT,at higher $\ell$). The spikes in TE are due to the denominator $C_\ell^{input}$ approaching to very small values.
import scipy as sp
from scipy import interpolate
from scipy import optimize
from scipy import integrate
from iminuit import Minuit
def get_delta_r(inputcl, lmax=200, wp = 7.6, fsky=1., n_steps=1024):
lmin=2
ell=pl.arange(lmin, lmax)
tenscl =interpolate.interp1d (ell,
pl.loadtxt('/Users/peppe/work/satellite_sims/Class_planck2015_r1_cl_unlensed.dat')[:,3][lmin:lmax])
lenscl =interpolate.interp1d (ell,
pl.loadtxt('/Users/peppe/work/satellite_sims/Class_planck2015_r0_cl_lensed.dat')[:,3][lmin:lmax ])
nl_w=interpolate.interp1d (ell, inputcl[lmin:lmax] )
clobs =lambda r,l :r *tenscl(l) + lenscl(l) + nl_w(l)
cltheor =lambda r,l: r *tenscl(l)+ lenscl(l) + nl_w(l)
logP =lambda r, r_inp,l: - 0.5*(2.*l +1 )*(cltheor(r_inp,l)/clobs(r,l) +pl.log(clobs(r,l)) -(2. *l -1)/(2.*l+1) *pl.log(cltheor(r_inp,l)))
#print "Looking for the maximum of L(r)--> "
func=lambda x: - pl.sum( logP(x, 0, ell) )
m=Minuit(func, x=1e-3 , error_x=0.0001, limit_x=(0.,1e-2), errordef=1e-8 , print_level=0)
m.migrad()
rmax= m.values['x']
L =lambda r: pl.exp ((logP(r,0,ell) ).sum() - ( logP(rmax , 0, ell) ).sum() )
Lmap= map(L, pl.linspace(0,.5e-2, 512) )/L(rmax)
a=0
b=1e-2
integr_steps=n_steps
Normalization= integrate.quad(L, a=a, b=b , limit=1000, epsrel=1e-8)[0]
R=0
dr= (b-a)/integr_steps
deltar=a
it_counter=0
while R<0.6801 and deltar<b :
deltar+=dr
R= integrate.quad(L, a=a , b=deltar, limit=1000, epsrel=1e-8)[0] / Normalization
it_counter+=1
deltar/=pl.sqrt(fsky)
#print "Delta r= %.2f x 10^-3 found after %d iterations at %g per cent C.L. "%( deltar*1000,it_counter, R )
return Lmap, deltar
l2=pl.arange(2, 501)
Lref,deltaref =get_delta_r( pl.ones_like(l2)*Nl(8.8) ,lmax=100 )
Luncorr,deltar1= get_delta_r( meanpre[2][2:lmax]*1e12 , lmax=100 , wp=wp )
Lcorr,deltar2= get_delta_r(meanpost[2]*1e12 ,lmax=100 , wp=wp )
fig=pl.figure(figsize=(10 , 5) )
pl.subplot(121)
pl.title('BB')
pl.loglog(l, cl[ 2:lmax+1,3] ,'k', label='lensing')
pl.loglog(l, meanpre[2][2:]*1e12 , label='pre-correction')
pl.fill_between(l, (meanpre[2][2:] -stdpre[2][2:])*1e12 , (meanpre[2][2:] + stdpre[2][2:])*1e12 , alpha=.4 )
pl.loglog(l, meanpost [2][2:]*1e12 , label='post-correction')
pl.fill_between(l,( meanpost[2][2:] -stdpost[2][2:] )*1e12 , (meanpost[2][2:] + stdpost[2][2:])*1e12 , alpha=.4 )
pl.loglog(l, pl.ones_like(l)*Nl0,'r', alpha=.5, linewidth=3,label='noise %.1f uK arcmin'%wp)
pl.ylim([1e-8,1e-4])
pl.xlabel(r'$\ell$', fontsize=16)
pl.ylabel(r'$\langle C_{\ell} \rangle_{MC} [ \mu K^2 ]$', fontsize=16)
pl.legend(loc='best')
foo=pl.xlim([2,600])
pl.subplot(122)
pl.plot(np.linspace(0,.5e-2, 512), Luncorr, label=r'$\delta r=%.2f \times 10^{-3} $'%(deltar1*1e3) )
pl.plot(np.linspace(0,.5e-2, 512), Lref, label=''+r'$\delta r= %.2f \times 10^{-3} $'%(deltaref*1e3),
color='red' ,alpha=.5, linewidth=3)
pl.plot(np.linspace(0,.5e-2, 512), Lcorr, label=r'$\delta r=%.2f \times 10^{-3} $'%(deltar2*1e3) )
pl.xlim([0., .5e-2 ])
pl.ylim([0.,1.1 ])
pl.legend(loc='best')
pl.xlabel(r'$r$', fontsize=15)
pl.ylabel(r'$\mathcal{L}(r)$', fontsize=15)
Caption: (left) BB power spectra from 20 MC sims, (solid) average power spectrum, (shaded area) the standard deviation. As a reference we plot as a thick red line the noise level of LB as expected after FG subtraction (assuming $7.6\, \mu K\, arcmin$) and the power spectrum of lensing B-modes. (right) Likelihood of $r$ to assess the level of contamination due to gain drift systematics in terms of $\delta r$ . Notice that the post correction likelihood (solid orange) is narrower than the reference one (solid red).